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Magnetic fields in the early universe could have played an important role in sourcing cosmological 
perturbations. While not the dominant source, even a small contribution might be traceable through 
its intrinsic non-Gaussianity. Here we calculate analytically the one, two and three point statistics of 
the magnetic stress energy resulting from tangled Gaussian fields, and confirm these with numerical 
realizations of the fields. We find significant non-Gaussianity, and importantly predict higher order 
moments that will appear between the scalar, vector and tensor parts of the stress energy (e.g., 
scalar-tensor-tensor moments). Such higher order cross correlations are a generic feature of non- 
linear theories and could prove to be an important probe of the early universe. 
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I. INTRODUCTION 
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It appears increasingly likely from cosmological data such as the cosmic microwave background (CMB) that the 
dominant source of perturbations in the universe was adiabatic and Gaussian in nature, consistent with the expecta- 
tions of an inflationary universe 0,0] • However, this does not exclude the possibility that non- linear sources might 
also have played a role in sourcing perturbations. While their impact on the power spectrum might be small, they 
| could dominate any non-Gaussianities that we observe. In order to search for these optimally we need clear predictions 
for the nature of the non-Gaussian signals that they would source. Here we begin to investigate the signatures for a 
' particular non-linear model, early-universe magnetic fields. 

Magnetic fields are observed on many scales in the cosmos, from planetary scales with fields of strengths of a few 
Gauss, to galactic scales at a strength of approximately /iG. They also exist between galaxies and are likely to be in 
clusters, both with field strengths somewhere between the nano- and micro-Gauss level and a coherence length on the 
1^ 1 order of megaparsecs. While fields on supercluster scales are extraordinarily difficult to detect, there are suggestions 
| that fields up to the order of [iG may exist even there. (See for example [E HI IE IE 0> HI f° r reviews.) 

The origin of these fields is uncertain. Possible creation mechanisms can be separated into processes occurring before, 
during or after recombination. Suggested post-recombination processes include the dynamo mechanism (I. : 9| (wherein 
magnetic energy is provided by rotational kinetic energy), or the adiabatic compression of an already-magnetized 
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cloud However, these still rely on a seed field, possibly created by some large battery mechanism; an efficient 
dynamo could amplify a field of strength 1CP 30 G to observed levels, while an adiabatic compression requires a seed 
many orders of magnitude larger. While mechanisms certainly exist to generate this field at reionization |l(J.lll| or at 
recombination itself 0,0], they might also have been created before recombination or even before nucleosynthesis; 
such early seed fields are constrained by limits from cosmic microwave background (see |I3 . ll5lllr^ |) and nucleosynthesis 

@iia 

Here we concentrate on fields from the early universe, and there is no shortage of suggestions for how such fields 
might have arisen. They could have been pro duced directly by inflation, or might have been generated during the 
electroweak symmetry breaking phase [E . 18, 19, 20. 21, 22, 23. 241. There are also recent studies into generic phase 
transitions generating large-scale fields (e.g.. |25ll26||). Cosmic defects might also be responsible [lEUE]- More recently, 
attention has been given to the possibility that the fields might have been created continuously in the period between 
lepton decoupling and recombination, through the vorticity naturally occurring at higher order in perturbation theory 
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Even after they are produced, they could evolve in the very early universe, perhaps as a result of 
hydromagnetic turbulence (e.g., |33j|). (Dolgov |34j provides a brief overview of many of these creation mechanisms.) 

Primordial magnetic fields can have a significant impact on the cosmic microwave background. While earl y tr eat- 
ments focussed on this effect for a Bianchi universe |35l l36| , more modern treatments [371 HE HE HE Ell l42l HE 
HI . HE HE HE HE HE HEl consider either small perturbations around a large-scale homogeneous field or a 'tangled' 
field configuration. Both of these scenarios source CMB perturbations, either directly through the scalar, vector and 
tensor stresses, or indirectly by the density and velocity perturbations they induce in the charged proton-electron 
fluids. Code for calculating the vector and tensor anisotropics generated by primordial magnetic fields was recently 
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added to the publicly-available CAMB [iM Hsij and that for scalars has been modelled independently poll4^ .l50 | . In 
addition to creating CMB anisotropics, magnetic fields can affect also CMB polarization by inducing Faraday rotation 

[EHESE1IE3. 

It would be useful to have a firm grasp on the statistical impact magnetic fields might have upon the microwave 
background. Our ultimate goal is to translate the non-Gaussianity in the magnetic stresses into a unique non-Gaussian 
signature in the microwave sky 56] . While there are many potential sources of CMB non-Gaussianity, including the 
non-linear evolution of the perturbations and gravitational lensing (for example, |5?1 l58t l59jl. there are also many 
ways a map may be non-Gaussian. The hope is that the different physical mechanisms will each have a characteristic 
non-Gaussian signature which can be searched for in the CMB sky. It is known that magnetic fields can introduce 
non-Gaussianity in many way s, including via Alfven turbulence or secondary perturbations, and each of these can be 
searched for in the CMB [60L l6ll l6^ | . In addition, helical magnetic fields can produce parity-odd correlations in the 
CMB polarization field [63. 

Here we focus on the non-Gaussianity of the magnetic sources themselves, with the three point moment as the start- 
ing point for our investigations. To this end we model the statistics of Gaussian-random magnetic fields, concentrating 
on the two- and three-point statistics of the magnetic stresses and the cross-correlations between the scalar, vector and 
tensor components, usually considered independently. We do so by both generating realizations of Gaussian-random 
magnetic fields and through numerically integrating pure analytical results, most of which are presented here for 
the first time. Significant non-Gaussianities and cross-correlations are expected due to the quadratic nature of the 
magnetic stress tensor; this suspicion is confirmed and the bispectra are presented. 

In section II we review the basic model of tangled primordial magnetic fields, their impact on the CMB and their 
statistics, detail our realizations of these fields and construct the stress-energy tensor. In section III we briefly consider 
the one-point statistics of the scalar pressures of the magnetic fields, evaluating the probability distribution function, 
the skewness and the kurtosis for both the iso trop ic and the anisotropic pressures. In section IV we review the two- 
point results (previously demonstrated in part |44|) and present the power spectra of the separate components of the 
stress-energy tensor, displaying the excellent agreement between the theory and the simulations. In section V we move 
onto the three-point moments, deriving and presenting analytical results for their bispectra and good agreement with 
the simulated fields. We discuss the implications of our results in section VI. 



II. TANGLED MAGNETIC FIELDS 



Rather than considering a large-scale homogeneous magnetic field with a small inhomogeneous perturbation, we 
instead consider entirely inhomogeneous fields tangled on some length scale (see for example Mack et. al. |44j). For 
simplicity (and to compare with the previous literature), we assume that the magnetic fields are random variables 
obeying a Gaussian probability distribution function; however, our realizations may be formulated more generally 
than this. Motivated by linear perturbation theory, we also take the fields to be effectively frozen, with their overall 
energy density decreasing as radiation (B oc a~ 2 , where a is the scale factor.) Due to the extremely high conductivity 
of the early universe [|| lisj , we also assume that the electric components of the electromagnetic field vanish. 

Magnetic fields, being a non-linear source with a full anisotropic stress, will naturally source scalar, vector and 
tensor perturbations. For scalar perturbations one does not expect a significant effect on large scales. However, the 
contribution to the temperature auto-correlation on the microwave sky can begin to dominate at an I of about 1, 000 
(see for example EH, H(J ) ■ This effect comes both from the impact of the magnetic energy and anisotropic pressure 
directly onto the spacetime geometry and from the Lorentz forces imparted onto the coupled proton-electron fluid. 
One also expects a significant impact on the vector perturbations as compared to the standard picture since the 
magnetic fields will both directly generate vector perturbations in the spacetime and will also contribute a Lorentz 
force to the solenoidal component of the velocity. Prior to neutrino decoupling a primordial magnetic field acts as 
a source for gravitational waves; following neutrino decoupling there will also be a contribution from the neutrinos, 
which Lewis has shown serves to cancel much of the magnetic stress |48j . Since the Lorentz forces and the stress-energy 
tensor are both quadratic in the magnetic field, we also expect a level of non-Gaussianity to be imprinted onto the 
fluid and perturbations, even for a magnetic field that is itself Gaussian. 

The features of the magnetic field - the Lorentz forces and the direct sourcing of geometric fluctuations - are 
all contained within the stress-energy tensor. We can then consider the statistics of the stress-energy tensor alone 
and be hopeful of characterizing the majority of the non-Gaussian effects that might impact on the CMB. The non- 
Gaussianity predicted from our analysis and our simulations can be projected onto the CMB sky by folding them with 
the transfer functions generated by the modified CAMB code |48|, |51| ; while this does not include support for scalar 
perturbations it will give good predictions on all scales for the -B-mode polarization correlations which is sourced 
purely by vector and tensor perturbations. However, care would have to be taken to detangle these from any £?-modes 
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caused by the gravitational lensing of the dominant -E-mode polarization. 1 



A. Statistics of the fields 

We begin by specifying the statistics of the tangled magnetic fields. In Fourier space, the divergence free condition 
for magnetic fields implies that 

(B a (k)B£(k')) - V(k)P ab (k)S(k - k') + ^H(k)e abc k c (2.1) 

where V(k) is the magnetic field power spectrum, P ab is the operator projecting vectors and tensors onto a plane 
orthogonal to k a and k bl 

P ab (k) = S ab - k a k b , (2.2) 



with 5 ab = diag(l, 1, 1) and TL{k) is the power spectrum of the anti-symmetric helical term (see for example 45, 63, 64].) 
Here we have assumed the fields are statistically isotropic and homogeneous. If the magnetic fields are Gaussianly 
distributed, then all their statistics are determined by their power spectrum. In the interests of simplicity we henceforth 
assume that the helical component of the field vanishes. 

In real space, the magnetic correlation function is the transform of the magnetic power spectrum, implying that 

(B o (0)B 6 (x)> = S ab C (x) + ——Ci(x) (2.3) 

where C (x) — J d 3 kV(k)e~ lkx and Ci(x) = j^js J d 3 k(V(k)/k 2 )e~ lk ' x . In the limit of very small separations, 
the correlation becomes 

(B o (0)B 6 (x)) = ^S ab C (0) + x a x b C 2 (0) (2.4) 

where C2(0) oc J dkk A V{k) is another correlation function related to the power spectrum. 
Often in the literature, the power spectrum is taken to be a simple power law, 

V(k)=Ak n . (2.5) 

To avoid divergences, these power law spectra are generally assumed to have some small scale cutoff associated with 
the photon viscosity damping scale; see Mack et. al. ^4| for an evaluation of these for vector and tensor perturbations 
generated by a stochastic field with a power-law spectrum. Durrer and Caprini 0, |6^ demonstrate that, for a 
causally-generated magnetic field, the spectral index must be n > 2; it must also in any event be n > —3 to avoid 
over-production of long-range coherent fields. They also show that, for the tensor case at least, the anisotropic stress 
will have a white-noise (n c g = 0) spectrum for all n > —3/2. They also demonstrate that nucleosynthesis limits on 
the gravitational waves produced by the magnetic fields place extremely strong bounds on mag netic fields, to the 
level of 10 _39 G for inflation-produced fields with n = 0, although this has been contested (|55, 66]). For spectra that 
might realistically imprint on the microwave background we must consider spectra with n < — 2, which are much less 
tightly constrained. For purposes of comparison we shall usually take either n = for a flat spectrum - where "flat" 
refers to the spectrum of primordial magnetic fields themselves - in which each mode contributes equally, or n = —2.5 
for a spectrum nearing the "realistically observable" n w — 3, which we shall refer to as "steep". Such fields can be 
produced by inflation (for example, fiH I22I] ) . 

It is worth emphasising that, while we are restricting ourselves to a power-law spectrum, this is not a necessity. 
We are generating simulated fields and calculating the statistics from these; it is as simple to employ other power 
spectra as it is to employ a power law. It is also a simple matter to employ a non-sharp damping scale - we could, for 
example, employ an exponentially or Gaussian-damped tail above a particular scale kp and gain some greater freedom 
in modelling the effective microphysics. As a particular example, Matarrese et. al. |30| derived the power spectrum 



1 As earlier noted, magnetic fields would also cause Faraday ro tati on from E modes into B modes lo2t l53t 154 | 55 l: support for scalar 
modes in a uniform magnetic field was added to CMBFast by \ r > j and also to CMBFast for tangled fields by l55l . both in the case of 
small rotations. 
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they expect from a field sourced from second-order vorticity in the electron-baryon plasma. This power spectrum has 
no simple functional form but is instead presented as a numerical integration and, in principle, it would be simple to 
employ this integration as the input power spectrum. 

The normalization of the power spectrum is typically fixed by reference to a particular co-moving smoothing scale 
A and the variance of the field strength at this scale, B\. Specifically, we smooth the field by convolving it with the 
Gaussian filter 

/(fc)=exp(-^J (2.6) 

and define the variance of the field strength at the scale A by 

B\ = (B a (x)B B (x)) , (2.7) 
implying that the power spectrum and B\ are related by 

B\ = J d 3 kV{k)e~ x2k2 . (2.8) 

This allows us to relate the astronomically observed field strengths at, say, cluster scales, to the amplitude of the 
magnetic power spectrum. 

For simplicity and concreteness, we will assume throughout that the fields themselves are Gaussian (consistent with 
previous literature). Obviously, any conclusions about the non-Gaussian signatures of magnetic fields will depend 
sensitively on this assumption, and we plan to explore more general scenarios in future work. Another ansatz for the 
fields is that they possess a % 2 probability distribution function; such fields might arise from weakly non-linear effects 
since they are sampled from the combination of two (first-order) Gaussian variables. For example, the sourcing of 
second-order vorticity in the electron-baryon plasma by first-order density perturbations could lead to a weak but 
continually generated magnetic field |29l l30L l3ll l32| . Magnetic fields sourced at recombination and later are unlikely 
to be reasonably described by the above formalism; those sourced at recombination |13| for example naturally arise 
at a very small scale. The question of the statistical nature of both recombination and very early universe fields has 
not been well explored in the literature. Here we present results from a Gaussian field as an illustrative example, but 
further work is needed to explore the consequences of any given microphysical mechanism. 



B. Realizations of magnetic fields 

To aid our study of the non-Gaussian properties of tangled magnetic fields, we create static realizations of the 
fields numerically. We create the fields on a grid in Fourier space of size Z| im , where i<jim is typically 100-200. The 
divergence free condition means we can generate the three magnetic field components for each k mode using two 
complex Gaussian uncorrelated random fields with unit variance, 

Ci 

c 2 

We then determine the magnetic field Fourier components by applying a rotation matrix, 

(B X \ 
B = I B y I = R C, 

where R is a 3 x 2 matrix. From the definition of the magnetic field power spectrum, we see that to get the proper 
statistical properties, we require 

(B a B b ) = R am (C m C n ) Rl b = (R • R T ) ah = P(k)P ab (k). 

While this does not specify the rotation matrix uniquely, it is straightforward to show that choosing the rotation 
matrix as 



R V(k) 1 / 2 



( k x k z ky 



k y k z ~ k x ] (2.9) 





FIG. 1: Sample realizations of the magnetic field for a spectrum n = 0. Left: Gaussian magnetic field slice, B x \ z =o; center: the 
(very non-Gaussian) isotropic pressure r| z= o; right: the (slightly non-Gaussian) anisotropic pressure r s | z= o. Compared to the 
magnetic field, non-linearity transfers power to smaller scales in the sources. 



will produce fields with the correct statistical properties. This rotation is well defined except in the case when 



k x = k y = 0. Here, B z — and the other components are uncorrelated, so we instead choose 



Rn 



V{k) 




(2.10) 



The reality of the fields is ensured by requiring B(— k) = B(k)*. 

Throughout, we are careful to avoid creating modes with frequencies higher than the Nyquist frequency of the grid, 
^Nyquist, which could be into power on other frequencies. Since the quantity of greatest interest, the stress-energy, is 
a quadratic function of the fields, it typically will have power up to twice the cutoff frequency of the magnetic fields. 
To avoid having aliasing of these fields, we generally require that the magnetic field cutoff frequency be less than half 
the Nyquist frequency. We also have an infra-red cut-off which is the inevitable result of working on a finite grid. 

Figure ^ shows a sample Gaussian realization of one component of the magnetic field along a slice through the 
realization, as well as the resulting trace and traceless components of the stress-energy (the isotropic and anisotropic 
pressures). Both the isotropic and anisotropic pressures show power on smaller scales than the fields themselves, 
a direct result of their non-linearity. In addition, the isotropic pressure is darker, reflecting a paucity of positive 
fluctuations and a significant deviation from Gaussianity. The anisotropic pressure appears to be more similar to the 
magnetic field - that is, relatively Gaussian. These observations will be made concrete in the next section when we 
consider the one-point statistics of the pressures. 



C. The stress-energy tensor 



Gravitational waves are sourced directly by the stress-energ y of a magnetic field, and the effect from the Lorentz 
force on the magnetized CMB also depends directly upon this |48| ; for this reason we shall concentrate our study on 
the space-space part of the magnetic stress-energy tensor, the stresses. In real space, these are 

r afc (x) = iBi(x)B((x)(f a6 - B a (x)B 6 (x). (2.11) 

In Fourier space, the product of two magnetic field components in real space becomes a convolution, 

f o6 (k) = J B a (q)B b (k~q)d 3 q. (2.12) 

The full stress-energy tensor in Fourier space is 

T ab (k) = ^S ab fu(k) - f ab (k). (2-13) 
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Even though the stress-energy is non-linear in the fields, we will be assuming throughout that the stress-energy 
perturbations in the magnetic fields are small compared to the total energy density. Thus, our treatment of the 
perturbations induced by the magnetic fields in the photons, baryons and so forth are purely linear. In this regime, 
the evolution of the scalar, vector and tensor perturbations decouple and we will thus decompose the magnetic field 
sources into these various pieces (interpreted as the isotropic pressure, anisotropic pressure, vortical and anisotropic 
stresses respectively) as 

Tab = ^ S "bT + {k a k b - 7^5 a b)T S 

+2k {a r^+rJ b . (2.14) 
We do this in Fourier space by applying combinations of projection operators; 

r(k) = <5 Qb T afc (k), 

r s (k) = (5 ah -(3/2)P ab (k))T ah (k)=Qab(k)T ab (k), (2.15) 
T <T( k ) = k (l P j)a {k)T %j {k) = V aii (k)ry, 

1 



r a J b (k) = I P a(i (k)P j)b (k) - -^(k)P ab (k) j TfcflO = T abij (k) Tij (k). 

In the above, (a ... b) denotes symmetrization in a and 6, i.e, A^ ab ) — (l/2)(A a b + Ab a ), and we are working with 
symmetrized projectors for the vector and tensor components for future ease. Here, hrf = hr? = t£ = 0. The first 
term in the stress-energy contributes purely to the scalar-trace part. For the others we can simply replace r a b with 

-Tab- 

The transfer functions which describe how fluctuations in the magnetic field stress e nergy im pact the microwave 
background have been previously evaluated, in various semi-analytical limits in [ifl ItH l42ll44j for the temperature 
perturbations, for example. However, the simplest route towards folding the expected non-Gaussianitics onto the 
CMB, will likely be from the CAMB code as modified recently by Lewis [i^,[^J; this produces the transfer functions 
generated by a magnetic field sourced before neutrino decoupling, neglecting the impact of the scalar perturbations. 
This could be supplemented by incorporating the extensions to scalar modes of Giovannini |49j into, for example, the 
CMBFast code We leave this issue to a later paper. 

III. ONE-POINT MOMENTS 

There are many ways to characterize non-Gaussianity, particularly given such a strongly non-linear stress-energy 
term. In this section we briefly consider the skewness and kurtosis of the one-point probability distributions of the 
isotropic and anisotropic pressures. In this section the results we present are the mean of twenty realizations with a 
grid-size of ld im = 192, and the errors quoted are one standard deviation. 

The simplest to consider is the distribution of the trace part, since it is simply the square of the magnetic field. 
Despite the divergence-free condition, the three components of the magnetic field at a single point in space are uncor- 
rected and Gaussian (as was shown above.) Thus we expect the trace of the stress-energy to have a \ 2 distribution 
with three degrees of freedom. 

All the moments of a one-point distribution may be given by its moment generating function as 

< = m = ^ i M(t)\ t=0 (3.16) 
where, for a x 2 distribution with p degrees of freedom 

M W = 7 ^72- ( 3 - 17 ) 

(1 - 2t) p/I 

The central moments are then readily calculated and the normalized skewness and kurtosis are defined to be 

= 72 = —2 -3 (3.18) 

We quickly find that, for the \ 2 distribution, the normalized skewness and kurtosis are 

7l = J- a 1.633, 72 = — = 4 (3.19) 
VP P 
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FIG. 2: The left figure shows the probability distribution of the isotropic and anisotropic pressures for a spectrum n — 0, with 
a Gaussian distribution shown for comparison. The isotropic distribution is well fit by a \ 2 curve. The right figure compares 
the anisotropic pressure distribution for different spectral indices. The x-axis is in units of the root mean square amplitude of 
the relevant field. 



where the numerical results are for a distribution with 3 degrees of freedom. The results from the realizations can 
be seen to be in agreement with the predictions; with a flat spectrum we find that, for the isotropic pressure, 
7i = 1.63 ± 0.01 and 72 = 3.99 ± 0.05. For a more realistically observable field, with a power spectrum of n = —2.5, 
say, we find 71 = 1.61 ± 0.01 and 72 = 3.92 ± 0.05. It is apparent that the statistics for the isotropic pressure are, as 
expected, relatively insensitive to the spectral index one employs. 

The anisotropic stress is harder to characterize because it is not a local function of the fields, but contains derivatives 
of them. However, it effectively is the sum of the products of two Gaussian fields which are, for the most part, 
independent of each other. The distribution of the product of two independent Gaussians is non-Gaussian but is 
symmetric (actually following a modified Bessel distribution, as shown in the appendix of Thus the effect of 

adding such terms is to dilute the skewness. That is, while the isotropic stress is the sum of three very skewed chi- 
squared distributed variables, the anisotropic stress is the sum of chi-squared terms and symmetric modified Bessels, 
making the result less skewed. 

The probability distributions of the isotropic and anisotropic stresses for a flat spectrum are plotted in the left 
panel of Figure [3 along with a Gaussian and a \ 2 distribution. The damping scale we employed was k c = /dim/4.1 - 
i.e. just beneath half the Nyquist frequency. For the steep spectra, we also used twenty realizations and a grid-spacing 
of Idim = 192 but with a damping scale at the size of the grid to ensure a reasonable mode coverage in the low- 
k/k c regime. The anisotropic pressure distribution has quite different properties when the spectral index is changed, 
including a switch in sign of the skewness. For the flat spectrum we find 71 = —0.24 ± 0.003 and 72 = 1.10 ± 0.01, 
while with a steep spectrum with a power spectrum of n = —2.5, we find 71 = 0.38 ± 0.01 and 72 = 0.86 ± 0.02. The 
distributions are plotted in the right-hand side of the figure, again with a sample Gaussian, and the change in the 
skewness is readily apparent. 

IV. TWO-POINT MOMENTS 

We next calculate the two point power spectra of the various perturbation types. These give a useful example of 
how the higher order calculations will proceed and give a means of testing our realizations. Some of these have been 
previously calculated, such as the vector [S[ and tensor [Til l4ll liif power spectra, while the trace and traceless scalar 
auto-correlations and their cross correlation, to our knowledge, have not. 

By the nature of the scalar-vector-tensor decomposition, we do not expect any cross correlations except between 
the trace and traceless scalar pieces. Thus, we consider five power spectra: one cross spectrum and the auto-spectra 
of the four pieces of the stress-energy. We focus on constructing rotationally invariant spectra which will contain all 
the information in the general correlations. 

In general, the power spectra will involve expectations of four magnetic fields. Since these are assumed to be 
Gaussian, they can be evaluated using Wick's theorem which, for four Gaussian fields, may be expressed as 



(ABCD) = (AB) (CD) + (AC) (BD) + (AD) (BC) 
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It is most useful to begin with the general two point correlation, 

<f a6 (k)f c * d (p)) = 5(k - p) J d 3 k'V(k')V(\k - k'l) (4.20) 

x (P ac (k')P bd (k - k') + P ad (k')P bc (k - k')) . 

(Note that there are two terms rather than three since we interested in the perturbations from the mean value of the 
field.) The power spectra of the various stresses may be obtained from this by applying different combinations of the 
projection operators (|2.15|) to yield, 

(r A (k)r B (p)) = S(k-p) J d 3 k'V(k')P(\k~k'\)T AB (4.21) 

with A and B denoting the two stress components and Tab = ^ab(j, Mi 0) denoting the relevant angular integrand. 
The relevant angles possible between the wavevectors have been defined as 

7 = k-k', fi = k'-k^k', /3 = k-k^k', (4.22) 

where k — k' denotes the unit vector in the direction of k — k'. 

The trace-trace correlation is found by applying the operator {l/4:)d ab d cd whence we obtain 

T Tr = \ (1 + M 2 ) . (4.23) 

Similarly, we obtain the traceless scalar auto-correlation function by applying (— l) 2 Q a fc(k)Q cc i(k); some algebra yields 
the result 

F T s T s = 2 + - | ( 7 2 + /3 2 ) - 3 7 M/3 + \l 2 P 2 - (4.24) 
The cross correlation between the trace and traceless scalar pieces requires the operator — (1/2)6 ab Q c d{k). This yields 

F TT s = -1 + | ( 7 2 + /3 2 ) + i M 2 - (4.25) 

For the vector and tensor contributions, it is useful to construct rotationally invariant combinations that can be 
relatively easily mapped onto the CMB. The divergenceless condition on the vectors implies that their correlation 
function can be written as 

(r^{k)rt v (p)> = \v v {k)P ab (\,)5(k p). (4.26) 

where our definition differs by a factor of two from Mack et. al. All the information is condensed in the rota- 
tionally invariant vector isotropic spectrum V v (k) — (r ( J / (k)T* 1/ (k)) . The operator necessary to recover this is 

k( a Pb)i(k)k^ c P^i(k) — k a k( c Pd)b{k) + kdk( a Pi>) c (k) . The resulting angular term can then be shown to be 

F T v T v = 1 - 2 7 2 /3 2 + m p. (4.27) 
Similar arguments apply for the tensor correlations. The full tensor two-point correlation is 

(rJ b (k)T c *J(p)) - \r T {k)M abcd {k)5{k - p), (4.28) 

where M a b c d{k) = P ac (k)P bd (k)+P ad (k)P bc (k) — P ab (k)P cd (k) which, as can be readily shown, satisfies the transverse- 
traceless condition on the tensors, and S ac 5 bd M abcd (k) — 4. We focus on the rotationally invariant tensor isotropic 
spectrum P T (k) — (r?(k)r* ) T (k)). Using the tensor projection operators and simplifying, the relevant operator 
is (P l{a (k)P b)l (k) - (l/2)P ij (k)P ab (k)) (P l(c (k)P d)l (k) - (l/2)P ij (k)P cd (k)) = P c(a (k)P b)d (k) - (l/2)P o6 (k)P cd (k). 
This leads to a simple angular term of 

F t t t t = (l + 7 2 )(l + /3 2 ). (4.29) 



The vector and tensor results differ from those otherwise presented (see Durrer et. al. and Mack et. al. |l7ll4ll l44| ). 
the vector case by 1 — j 2 as a result of employing the symmetrized projection, and the tensor case by 7 2 — (3 2 . 



9 



8 



6 



2 













o 


<TT> 








S S 






□ 


<T X > 








S 

<TT > 
















< T aV> 


- 






<x ab Tx ab T> 




^^t^-tr*^-^. | .T^T 1 !. 









1 2 

kkc 

FIG. 3: Magnetic stress-energy power spectra (n = 0) - the realizations agree well with the analytic predictions. The error bars 
from the realisations are small except at very low k. 



It is straightforward, however, to see that if one redefines the integration wavemode as k' = k" — k one maps 
fx' — > /i, j3' — y 7, 7' — y f3. On integration, then, the product 7 2 /3 2 is invariant while 7 2 — (3 2 may be taken to vanish. 
Our results are thus in agreement with those previously presented. 

We can compare numerical integrations of these power spectra with the results arising from the realised magnetic 
fields. Our results for a flat power spectrum (n = 0) are presented in Figure |3 where P(k) denotes the various spectra 
(all presented on the same scale). The agreement between the analytic results (lines) and a simulated field (data 
points) is striking. We have plotted the power spectra averaged over twenty realizations with a grid-size of /dim = 192, 
a damping scale at k c = i<ji m /4.1, with error bars of l-a. We have also rebinned the results because otherwise the data 
points obscure the theory. 

We could also consider magnetic fields with n > 2 corresponding to causally-generated fields ^?|k the features for 
such fields are qualitatively similar to those for a flat spectra and there are no difficulties in evaluating the theoretical 
predictions in this regime. However, as commented, Caprini and Durrer demonstrated that primordial fields of this 
type would be unobservably small in order to not violate nucleosynthesis bounds on gravitational waves. Moreover, 
blue spectra naturally tend to pile power around the cut-off scale k c and since we are not at all modelling the 
microphysics the results would have to be treated with caution. 

Analytic results can be found in certain limits. There are two regimes of interest for the spectral index, as shown 
by Durrer et. al. |4l|]. For n > —3/2 the integrations are dominated by the cutoff scale, resulting in constant spectra. 
In this regime, if k <C k c , the angular integrations are straight forward (fi ~ —1,(3 ~ —7.) Relative to the trace 
correlation, the amplitudes of the other correlations are t s /t — 7/5, r x /r = 0, r v /r = 14/15 and t t /t = 28/15 
respectively. For n < —3/2 the situation is considerably more complex and we content ourselves with the results of 
our simulations in Figure^ Immediately apparent is that the effects of the cutoff are reduced by the strongly tilted 
magnetic spectrum, with each spectrum quickly approaching the power law, Pa(&) oc k 2n+3 , that is naively expected 
from the fc-integration. Also notable is the change of behaviour of the scalar cross-correlation; whereas this vanishes 
on large scales in the n > —3/2 regime, it remains finite on large scales for n < —3/2 and so in principle might be 
observable on the sky. There are also the effects of the infrared cut-off causing a suppression at low-fc. 

V. THREE-POINT MOMENTS 

In this section we focus on the three point moments in Fourier space, for which it is possible (if laborious) to obtain 
analytic expressions. The magnetic field realizations provide a way of exploring other kinds of non-Gaussianities which 
may arise. 

At the three point level, it is no longer guaranteed that correlations between the scalar, vector and tensor pieces 
will vanish, and we present some of the first calculations of these here. There are many possible three point moments, 
but here we consider only the rotationally invariant combinations (ttt), (ttt s ), (tt s t s ), (t s t s t s ), {tt^t^) and 
(t t„ t„ ). Due to their complexity we postpone calculations of the correlations with the tensor components to a later 
paper and content ourselves with presenting the results from the realizations; for completeness these are (ttJ^t^) and 
{T S T^ b T^ b ) for correlations with the scalars, {Ta T Jb T b) the vectors and the auto-correlation (t^t^t^). We work 
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FIG. 4: Magnetic stress-energy power spectra (n = —2.5). The dashed line shows the spectral dependence expected naively, 
oc k 2n+s . The turnover at low k reflects the finite size of the grid. 



throughout in Fourier space, where the three-point moments are known as the bispectra. One advantage of working in 
Fourier space is that the transfer functions, which fold in the fluid dynamics and describe the impact on the microwave 
background, are local. 



A. General considerations 



In principle we can calculate all the three point statistics described above in Fourier space. The bispectra involve 
three wave modes, and since we assume the fields are homogeneous and isotropic, the sum of the three modes must 
be zero. Thus the bispectra are a function of the amplitudes of the modes alone (or alternatively, two amplitudes and 
the angle between them.) We denote different geometries by selecting a baseline k and a vector p making an angle <f> 
with k and having an amplitude p = rk (see Figure EJ. We may then calculate q = — k — p. For simplicity, we here 
concentrate on the colinear (degenerate) case in which p = k implying q = —2k - that is, r = 1 and = 0. 

We calculate the bispectra analogously to the power spectra, although matters are complicated by the need to 
deal with expectations of six fields rather than four. The object of most general interest is B%j uimn (k, p, q) = 
{fij (k)ffe;(q)r mn (p)) , which is related to the expectation value of six magnetic fields, 

B ljklmn (k,p,q) = JJJ dWp'd 3 q' (B i (k')B i (k - k')P fc (p')P;(p - p')B m (q')B n (q - q')> . (5.30) 

As in the two-point case, all three-point moments of interest may be found by applying the relevant projection 
operator, Aijkimn, to this. Expanding this six-point correlation with Wick's theorem generates fifteen terms, eight of 
which contribute to the reduced bispectrum, that is, the bispectrum neglecting the one-point terms proportional to 
<5(k), <5(p) or S(q). This leads eventually to 

Syfeimn(k,p,q) = (5.31) 
<5(k + p + q) J d 3 k'V(k')T(\k - k'|)7>(|p + k'|)P ifc (k') (P im (k - k')P„(P + k')) + {* ^ j,p -> q} , {k ~ I, m <-> n.} 

These eight terms reduce to the same contribution if the projection tensor Aijkimn that recovers a set bispectrum is 
independently symmetric in {ij}, {kl} and {mn}. 

In the power spectra calculations, the geometry was straight forward; here it is considerably more complicated. 
The three wavevectors of the bispectrum are constrained by homogeneity to obey k + q + p = 0. Combined with 
the dummy integration wavevector, these define a four sided tetrahedron. This has six edges, k, q, p,k',k — k' and 
p + k'. From these, we can generate fifteen unique angles which could come into the bispectra calculation. This is 
to be compared to just three edges and three angles required for the power spectra. Clearly these angles are not all 
independent; they are, in fact, functions of just five underlying angles. 

For our purposes, it is easiest to work with the fifteen which we separate into four hierarchies; those between the 
set wavevectors k, p and q, angle cosines of these vectors with k', cosines with k — k', and cosines with p + k'. The 




final group are defined below. We take the angles, with {a, b} C {k, p, q}, to be 

ob = a-b, a a = a-k', /3 =ak-k', 7 Q =ap + k', 

(3 = k' ■ k^k' 7 = k' ■ pTk' p = k^k' ■ p + k'. (5.32) 

In terms of the angles ^k q and ^ pq in Figure [S] we obviously have 9k q = — cos^ 9 and similarly for pq . We also have 
afc = cos 8. 

In a manner entirely analogous to the two-point case we find the different bispectra by applying to l|5.31[l different 
projection operators to extract the relevant scalar, vector or tensor parts and, given that we ensure that Aijkimn has 
the required symmetries, we may express the bispectra as 

(TA(k)rs(q)Tc(p)) =£(k + p + q) j d 3 k'V(k')V(\k - k'\)V{\p + k'|) {SFabc) (5.33) 

where {ABC} denote denote different parts of the stress-energy tensor and Tabc is the relevant angular component. 

B. Scalar bispectra 

We begin with the simplest case, the bispectrum of the magnetic pressure. This is found by defining Aijkimn — 
(l/8)6ij5ki5 mn which gives us 

8F TTT = ^ + f + tf- pjp. (5.34) 

The first scalar cross-correlation will be between the square of the pressure and the anisotropic pressure, found by 
using Aijkimn = {-l/4)8ij5kiQ m n(<l) to give 

-8T TTT s = 3 (l - a\ - - 0* - i(/5 2 + f + p 2 ) + a q {P q ~P + l q l) + fi{P q l q + ~ PlM^j (5-35) 

Similarly the second scalar cross-correlation, with Aijkimn — 0-/2)SijQkl(p)Qmn( ( l)i gi yes 

5 

T TT s T s = Y, :F ?tStS ( 5 - 36 ) 

n=0 
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where 

8^ r r s rS = -6, 
8J 7 ^ T s T s = 0, 

8T 2 TT s T s = W + f + fi 2 + 3 (a 2 p + a 2 + ft + (i 2 q + 7 2 + 7g 2 ) + 9^, (5.37) 
8J> T s T s = - (/?7M + 3/2(/3p7 P + /3,7,) + 7(a P 7 P + "97?) + /3(«p/3 P + £*<A) + 99 pq (a p a q + /3 p f3 q + j P i q )) 
8J^ tStS = 3 (P(p,a p j p + 7/Vfg + 3a p (3 q 6 pq ) + 3 (a P 7 P a,7, + /3 P 7 P /3 g 7 g )) 
8^ r s r s = -9f3a P i P f3 q j q . 

Finally the anisotropic scalar bispectrum is found by applying Aij km i n = (—l) 3 Qij(k)Q k i(p)Q mn (q) which results 

in 

6 

T T s T s T s = ^ !F™s T s T s (5.38) 

71=0 

with 

— 8J-®s T s T s = 9 
— 8J-^.s T s T s = 

-8^ t 2 StStS - - (^ 2 + 7 2 + f + 3(a 2 +a 2 p + a 2 q +p 2 +Pl+ f3 2 + 7fe 2 + 7 P + 7,) + 9(^ p + 9 2 q + 9 2 pq )) 
-8^s r s T s = 3 ^A(/3fe7fe + /?p7p + A?7<? + ^7) + 7(afc7fc + a p -f p + a q j q ) + (3(a k (3 k + a p (3 p + a q (3 q ) 

+39 kp (a k a p + f3 k f3 p + j k j p ) + W kq (a k a q + f3 k f] q + j k j q ) + 36 pq (a p a q + P p f3 q + 7^) + 99 kp 9 kq 9 pq 

-8T^s T s T s = -3 (jp,a k f3 k + (3fia p j p + (3jf3 q j q + 3(a k f3 k (a p f3 p + a q (3 q ) + a p j p (a k ~/ k + a q j q ) + f3 q "/ q ((3 k -f k + f3 p j p )) 

+3(/l6> fc p/? fc 7p + ^9 kq a k ^ q + (39 pq a p {3 q ) + 9(9 kp 9 kq j p j q + 6 k p6 pq (3 k (3 q + 9 kq 9 pq a k a p ) 



-8J*s T s T s = 9 (^a k p k a p j p + ja k f3 k f3 q j q + (3a p j p {3 q -f q + S(9 kp f3 k ^ p (3 q j q + 9 kq a k a p j p j q + 9 pq a k f3 k a p /3 q )) 
-8J^s T s T s = -27a k f3 k a P 7 P f3 q -f q . 

C. Cross bispectra 

For the vector and tensor correlations we restrict ourselves to the various rotationally-invariant quantities, which 
can be identified with cross-correlations between the scalar pressures and either the vector or tensor moduli. The 
first of these, the correlation between the scalar pressure and the vorticity, we recover with the operator Aij k i mn — 
(l/2)5ijP( fe Pj) (p)g( m P n ) (q). The eventual result is 

6 

T TT v T v = Y J ^ r y (5-39) 



with 



71 = 1 



8J- TT v T v — —39 pq 

8T 2 T v r V = Ip^q 

SJ* T v T v = M(/3 P 7 9 + Pqlp) + 7(a P 7g + a 9 7 P ) + 29 pq (a 2 p + (3 2 + 7 2 + a 2 + (3 2 + 7 2 + i/3 2 + 29 2 pq ) 

8f* T v r v = ~P{fia plq + 77 P /3 g ) - 206 pq (a q p q + a p [3 p ) - 2 lplq {a 2 p + (i 2 p +a 2 q + (3 2 + \fi 2 + 29 2 pq ) 

-2(a p a q + I3 p f3 q ){-f 2 + 7 2 + 29 2 pq ) 

8J* T v T v = 49 pqlplq (a p a q + (3 p (3 q ) + 2{3 lplq {a p f3 p + a q (3 q ) + 2^a p /3,( 7 2 + 7 2 + 29 2 pq ) 

8J^ T v T v = -4:/39 pq a p -f p (3 q j q . 
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The cross-correlation with the anisotropic pressure is recovered with the operator Aijkimn = 
a(p)Q(m-Pn)a{ ( l)) the eventual result is 

7 

T T s T v T v =-J2KsrV T v (5-40) 

n=l 

with 

8T\s r v T v — 69 pq 

-8^s T v T v = 4 7p 7 9 

-8^s T v T y = (PpJ q + /3 q Jp) ~p + {a p ~f q + a g7p ) 7 + (39 kp J q + 39 kq J p ) 7fc 

+9 pq (69 2 kp + 69 2 q + A9 2 pq + 3a\ + 3/3 2 + 2« 2 + 2(3 2 + 2 7 2 + 2a 2 q + 2(3 2 + 2 7 2 + (3 2 ) 

8J^s T v T v = A9 pq {39 kp 9 kq + CtpCtq + Ppfiq + 7p7l? ) + f39 pq (3Q!fc/3 fe + 2 (otp(3p + C(q(3 q )) 

+69 pq {Okp (a k a p + flkflp) + 9 kq (a k a q + f3 k f3 q )) 

+l P lq (3a 2 + 3(3 2 + 2a 2 p + 2/3 2 p + 2a 2 q + 2(3 q + 69% + 69 2 q + f) 

+ 2 (7p + lq) (a p a q + (3 p (3 q + 39 kp 9 kq ) + f3 {a p ^ q Ji + /3«j 7 p 7 ) + 3 {a k j p ^9 kq + /3 k j q Jl9 kp ) 

+3 7/c (a k a p j q + f3 k P q j p ) 
-8T^s T v T v = 49 pq (a p /3 q P + 39 kq a k a p + 39 kp f3 k p q ) + 4 (9 pq (a p a q + (i p (3 q + 39 kp kq ) 7p7 

+3a k f3 k (a p (3 p + a q (i q )) + 69 kp ((a k a p + (3 k j3 p ) 7p7? + (3 k f3 q ( 7 2 + 7 2 )) 

+69 kq ((a k a q + /3 k f3 q ) 7p79 + a k a p ( 7 2 + 7 2 )) + (la p /3 p + la q (3 q + 3a k (i k ) j p j q /3 

+2a p /3 q (3 ( 7 y + 7 2 ) + 3a k f3 k (apjqji + A? 7f)7 ) 
8T®s T v T v = 129 2 q a k apl3 k p q + A9 pq j p j q (a p P q P + 39 kq a k a p + 39 kp /3 k P q ) 

+6a k f3 k ~f p ~f q (a p /3 p + a q (3 q ) + 6a k f3 k a p f3 q ( 7 2 + 7 2 ) 
-8T T s r v r v = I29 pq a k a p (3 k p q j p ^ q . 

The full tensor correlation, (t^t^t^), has not been calculated, but it can be found by the application of 

AijUmn = (P o(i (k)P j)6 (k) - I^(k)P a fc(k)) (P b ( fc (p)P)c(p) - iP fei (p)Ac(p)) (Pa(m(q)P„)c(q) - ^P m n(q)Pa C (q)) . 

If we consider first the colinear case, this reduces immediately to a product of projection tensors on k. Employing 
the symmetries of Aijkimn in {ij}, {kl} and {mn} and of Bijkimn m {ik}, {jm} and {In} one may then demonstrate 
that this vanishes identically. The general case is not straight forward and we defer it to a later study. 

D. Flat spectrum results 

As was the case for the power spectra, there are two very different spectral regimes for the bispectra. For flat spectra, 
n > —1, the integrals are dominated by the highest k modes around the cutoff scale. For these spectral indices, the 
bispectra become independent of k when k ^ k c and the analytic expressions are straight forward to integrate. 
Indeed, we can do them exactly in the limit k <C fc c ; we present these results below both for a generic geometry and 
for the colinear case. For convenience we define the generic geometry using the two angle cosines a k — cos((9) and 
9 kq = — cos(^. g ) from which our specification of r and (\> may be recovered - see Figure 

In the general case, setting B = A 2 k^- n+1 ^ /3(n + 1) with A the amplitude of the magnetic field power spectrum, 




we find 

(r(k)r(p)T(q)> = BwS(k + p + q) 
(r(k)r(p)r s (q)) = 

(T(k)r 5 (p)T S (q)) = ^Btt (2 + 6cos(0)cos(&,) - 3cos 2 (a 9 ) - 3 cos 2 (0) + 6 cos 2 (0) cos 2 (&,) 
-6cos(</>) cos 3 ^) - 6cos 3 0) cos(£ kq ) + 6cos 3 (</>) cos 3 (^ 9 )) S(k + p + q) 
(r s (k)T S (p)r s (q)) = (l - 3 cos 2 (0) cos 2 ^) - 3 cos((/») cos(^ 9 ) sin 2 (0) sin 2 (6 9 )) <5(k + p + q) (5.41) 

( r ( k ) T ( T(p) T ( r(q)> = -^ BlT (3cos(0)cos(^ fc9 ) +sin 2 ((/))sin 2 (^ 9 ) - 3 cos 3 (0) cos(£ fcg ) ~ 3 cos(0) cos 3 (£fc g ) 

+4cos 3 (0) cos 3 (^ fe g) - sin 2 (0) cos 2 (0) sin 2 (5^) - sin 2 (</>) sin 2 (^ fe(? ) cos 2 (^ 9 ) 

+4 cos 2 (0) cos 2 (£ fc9 ) sin 2 (</>) sin 2 S(k + p + q) 
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( T ' S ( k )' r i r(p)^(q)) = Yo5 B?r ( 6cos ('? !, ) co s(Cfeg) - 6cos 3 0)cos(£fc 9 ) - 6cos(0)cos 3 (Cfe 9 ) - sin 2 (0)sin 2 (£ fe9 ) 
+8cos 3 (0) cos 3 (£ fc(Z ) - 2sin 2 (</>) sin 2 (£ fe(Z ) cos 2 ((, kq ) - 2sin 2 (0) sin 2 (£ feg ) cos 2 (0) 
+8cos 2 (0) cos 2 (a g ) sin 2 ((/>) sin 2 (a,)) <5(k + p + q) 

Specialising these to the colinear case wherein (f> = S, k q — these reduce to 

(r(k)r(p)r(q)) = BwSQt + p + q) 
(r(k)r(p)r s (q)) = 

(r(k)T S (p)r S (q)) = I^(k + p + q) 

(r s (k)r s (p)r s (q)) = -^^(k + p + q) (5.42) 

(r(k)r v (p)r v (q)> = -^^(k + p + q) 
^4 

(r^k^p^q)) = — i?7rJ(k + p + q) 

The bispectra that are derived from the simulated fields are heavily compromised by the grid-size; unlike the two- 
point case the three-point moments use only a restricted number of the modes, selected by the geometry chosen for 
the wavevectors. The result from a single realization is in most cases noise-dominated. To overcome this difficulty 
we have chosen to simulate a large number of different realizations, taking the mean signal and using their variance 
to provide an estimate for the errors involved. The results, for a flat power spectrum, a grid-size of Zdim = 192 
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FIG. 7: Here we show the colinear bispectra for a steep spectrum (n = —2.5.) The scaling is as expected naively, oc fc 3n+3 = fc 9 / 2 
(dotted line). 



(and a damping scale of k c = Zdim/4.1) and 1, 500 combined realizations, are plotted rebinned in Figures El with the 
numerically-integrated predictions overlaid. For simplicity we have concentrated on the colinear case, wherein p = k 
and so q = — (k + p). We plot k/k c against B{k) where B(k) represents the colinear bispectra. 

E. Steep spectrum results 

In the regime n < — 1 matters are, as with the two-point case, complicated by the presence of numerous poles; 
the integrals are dominated by the volume lying between the poles. Rather than attempt a solution, we use our 
realizations to calculate the bispectra. We present the results for n = —2.5 in Figures [7] There is a dependence on 
the gridsize for this spectrum due to the paucity of modes of appreciable power given the strongly red spectrum. 
We tested this by running three realizations with differing grid-sizes, keeping the total number of modes constant in 
each case constant; specifically we took 1,500 simulations at /dim = 192, 5000 simulations at Idim = 128, and 40,000 
simulations at Idim — 64. As might be expected we found a suppression for the case wherein Idim = 64 - due to 
scarcity of modes at low k - but there was good agreement between the other two cases. We have plotted the results 
from the Idim — 192 run for greatest dynamic range, again rebinning into 64 bins. 

With this highly-tilted spectrum we see, as with the two-point case, that the features of the magnetic spectrum at 
the cut-off scale apparent in the flat case are washed out by the spectral tilt. The predominant shape again is the 
fc-dependence we expect from the radial component of the integral; in this case B(k) oc fc 3 ( n + 1 ) . (These are plotted for 
comparison.) Again, as with the two-point case there is an infrared suppression. The magnitudes of the bispectra fall 
into three close bands; the strongest are the correlations between the scalars and the tensors, while the middle-band 
is composed of the scalar auto- and cross-correlations. The weakest bispectra are those involving correlations with 
the vectors and, in some cases, the 1 — a error bars are greater than the mean value; such points have been removed 
from these plots for aesthetic purposes. The (t^t^t^) correlation is particularly weak. 



VI. DISCUSSION AND CONCLUSIONS 



We have studied, analytically and via realizations, some of the higher point correlations for tangled large-scale 
magnetic fields. The analysis is obviously highly model-dependent, and extending the analytic results to other models 
could be very difficult. However, it would be a simple matter to generalize the numerical realizations to perform 
the same analysis for a wide variety of time-independent models; all we require is the statistical distribution of the 
underlying field, the power-spectrum and the form of the stress-energy tensor. For example, while we have assumed 
the simplest - and perhaps most instructive - case of a magnetic field with Gaussian statistics, it would be a simple 
matter to employ a \ 2 probability distribution for a magnetic field, which is physically motivated by the creation 
mechanism considered by Matarrese et. al. |30j |. (Since our code is time-independent the interpretation would be of 
the final field immediately prior to recombination.) The realizations will however be limited by the narrow dynamic 
range allowed in the computation. 
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In our particular case of a tangled primordial magnetic field with a power-law spectrum, we have demonstrated 
that we can recover the 1-, 2- and 3-point statistics from simulations and with an excellent agreement to our analytic 
predictions. At the one-point level we have not only verified that the magnetic energy density follows a \ 2 probability 
distribution function (as expected given that it is directly the square of the underlying Gaussian magnetic field), but 
that the anisotropic pressure is also non-Gaussian with a significantly more complex relation to the fields. There is 
also a spectral dependence on its probability distribution function affecting the skewness, which swaps sign as one 
passes through n = —3/2. The kurtosis, while exhibiting a spectral dependence, remains positive. 

At the two-point level we have calculated the scalar, vector and tensor auto-correlations, as well as the scalar-cross 
correlation. For the scale-invariant spectrum we confirm the power-spectra with the expected ratios for scales longer 
than the cutoff; we also see that the scalar cross-correlation vanishes on large scales but is not in general entirely zero 
for modes close to the cutoff scale. For a highly-tilted spectrum, the cutoff scale is less important and the spectra 
behave with power law behavior and with the relative ratios approximately constant. The fc-dependence is the power- 
law that would be expected from a naive point of view. The surprise is that in this regime the scal ar cross-corr elation 
no longer vanishes on large-scales; rather, the correlation remains roughly constant at (rr s ) / ' \J {tt){t s t s ) w 0.7. 
Thus, while the isotropic and anisotropic pressures are indeed correlated, they are not perfectly correlated. 

At the three-point level we considered a number of rotationally-invariant bispectra, concentrating for simplicity on 
the colinear case. We find significant non-Gaussianities - in excellent agreement with the analytic predictions. For the 
flat spectrum, these approach fixed ratios on large scales. These can be both positive or negative, or even zero. As 
in the two-point case, some qualitative aspects change when we consider a strongly-tilted magnetic power spectrum; 
features arising from the cutoff scale tend to disappear, leaving instead a simple power-law drop off. The relative 
ratios of the bispectra also change, and even their relative signs differ. The correlation between the isotropic pressure 
squared and the anisotropic pressure, which vanishes for scale invariant spectra, becomes non-zero. 

It remains for future work .56] to fold these results in with the transfer functions for magnetized cosmologies and 
calculate the non-Gaussianities expected to be imprinted upon the cosmic microwave sky. We can then consider how 
such signals might be used to constrain the properties of a magnetic field of this type. It would also obviously be 
straight-forward to consider different magnetic power spectra and statistics. While it is too early to speculate what 
these will discover, the higher order correlations in the sources will certainly lead to similar higher order correlations 
in the CMB observables, including perhaps cross correlations between the polarization modes such as (E 2 B 2 ). 

The techniques used in this paper could be applied to a broad variety of models. Moreover, we have here to 
our knowledge presented the first calculations of cross-correlations between, for example, scalars and tensors, and 
demonstrated that they can be of an equal magnitude to scalar auto-correlations. This has great potential relevance 
to the wider field of sources with non-linear stress-energy tensors, or sources with non-Gaussian initial conditions. For 
example, it would be interesting to compare our results to the same moments evaluated for defect models, particularly 
given the current resurgence of interest into networks of cosmic strings. 
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